Code
library(here)
source(here("02_scripts/03b_fun_calib_ntd.R"))
theme_set(theme_minimal())
opt_params_notimedepend <- readRDS(here("01_data/map_calib_ntd.RDS"))Comparison between MAP calibrated model outcomes and calibration targets - not time dependent
Used overall deaths, prevalence of rx opioid use, and total OAT as calibration targets
Code
mod_opt_map_ntd <- mod_calib_ntd(as.numeric(opt_params_notimedepend$par))1 Prevalence of Rx opioids use
Code
# prevalence of people on prescribed opioids
prop_opioids_rx_calib <- prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = length(2015:2018),
list(NULL,
paste0("prop_opioids_rx_",
str_sub(2015:2018, 3, 4)))))
tot_pop_calib_tbl <- tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
nrow = 12,
ncol = length(2015:2018),
list(NULL,
paste0("tot_pop_",
str_sub(2015:2018, 3, 4)))))
for (i in 1:length(2015:2018)) {
yr <- 2014 + i
prop_opioids_rx_uncalib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_uncalib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
prop_opioids_rx_calib[, i] <- assign(
paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
"BPO_OTHER", "BPO_PALLIATIVE")]) /
rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
tot_pop_calib_tbl[, i] <- assign(
paste0("tot_pop_", str_sub(yr, 3, 4)),
rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
)
}
prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(2015:2018,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_uncalib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_calib_tbl <- prop_opioids_rx_calib %>%
pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>%
arrange(grp) %>%
mutate(year = rep(2015:2018,
each = 12),
year_mon = paste(year, month.abb[1:12], sep = "_")) %>%
bind_cols(., tot_pop_calib_tbl %>%
pivot_longer(1:4, names_to = "tot_pop",
values_to = "tot_pop_val") %>%
arrange(tot_pop) %>% select(-tot_pop))
prop_opioids_rx_target_tbl <- calib_target_tbl %>%
filter(group == "prev_on_opioidsrx") %>%
select(year, target, group) %>%
rename(grp = group,
target_val = target) %>%
mutate(year_mon = paste(year, month.abb[6], sep = "_"))
prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Uncalibrated Model") %>%
bind_rows(., prop_opioids_rx_target_tbl %>%
select(year, target_val) %>%
mutate(grp = "target")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(grp = "Calibrated Model"))
prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_calib_tbl$year_mon <- factor(prop_opioids_rx_calib_tbl$year_mon,
levels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"),
labels = paste(rep(2015:2018,
each = 12),
month.abb[1:12],
sep = "_"))
prop_opioids_rx_tbl <- prop_opioids_rx_uncalib_tbl %>%
mutate(group = "Uncalibrated Model") %>%
bind_rows(.,prop_opioids_rx_calib_tbl %>%
mutate(group = "Calibrated Model"))
wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>%
mutate(grp = "Target") %>%
bind_rows(., prop_opioids_rx_uncalib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
grp = "Uncalibrated Model")) %>%
bind_rows(., prop_opioids_rx_calib_tbl %>%
group_by(year) %>%
summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>%
ungroup() %>%
mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
grp = "Calibrated Model"))
prop_opioids_rx_tbl$group <- factor(prop_opioids_rx_tbl$group,
labels = c("Calibrated Model", "Uncalibrated Model"),
levels = c("Calibrated Model", "Uncalibrated Model"))
p1 <- ggplot() +
geom_point(data = prop_opioids_rx_tbl, aes(x = year_mon, y = target_val,
color = group, fill = group)) +
geom_point(data = wprop_opioids_rx_tbl %>% filter(grp == "Target"),
aes(x = year_mon, y = target_val, color = grp, fill = grp), size = 2) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
xlab("Year_Month") + ylab("Prevalence of Rx opioid use") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
plotly::ggplotly(p1)1.0.1 Comparison between uncalibrated annual weighted average and observed targets of prevalence of Rx opioids use
Code
wprop_opioids_rx_tbl$grp <- factor(wprop_opioids_rx_tbl$grp,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p2 <- ggplot() +
geom_point(data = wprop_opioids_rx_tbl,
aes(x = year_mon, y = target_val,
color = grp, fill = factor(grp)), size = 2) +
xlab("Year") + ylab("Prevalence of Rx opioid use") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# ylim(0,0.2) +
plotly::ggplotly(p2)2 Deaths and OD-Deaths
Code
################ Deaths
# number of deaths
num_deaths_calib <- num_deaths_uncalib <- rep(NA, length(2016:2020))
for (i in 1:length(2016:2020)){
yr <- 2015 + i
num_deaths_uncalib[i] <- mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
num_deaths_calib[i] <- mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"] -
mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_DEATH"]
}
# number of OD-deaths
num_od_deaths_calib <- num_od_deaths_uncalib <- rep(NA, length(2016:2021))
for (i in 1:length(2016:2021)){
yr <- 2015 + i
num_od_deaths_uncalib[i] <- mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
num_od_deaths_calib[i] <- mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"] -
mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
year_mon_cycle_tbl$mon == 12] + 1,
"BO_OD_DEATH"]
}
deaths_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_deaths", "total_od_deaths")) %>%
mutate(target = ifelse(target == "total_deaths", "Total deaths",
ifelse(target == "total_od_deaths",
"Total OD-deaths", NA)),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_deaths_uncalib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
mutate(target = "Total OD-deaths")) %>%
mutate(year = c(2016:2020, 2016:2021),
group = "Uncalibrated Model")) %>%
bind_rows(., data.frame(target_val = num_deaths_calib) %>%
mutate(target = "Total deaths") %>%
bind_rows(., data.frame(target_val = num_od_deaths_calib) %>%
mutate(target = "Total OD-deaths")) %>%
mutate(year = c(2016:2020, 2016:2021),
group = "Calibrated Model"))
deaths_target_uncalib_tbl$group <- factor(deaths_target_uncalib_tbl$group,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p3 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# +
# facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(p3)Code
p4 <- ggplot() +
geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total OD-deaths"),
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total Opioid-related Overdose Deaths") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
plotly::ggplotly(p4)3 OAT
Code
# number of oat
num_oat_calib <- num_oat_uncalib <- rep(NA, length(2018:2021))
for (i in 1:length(2018:2021)){
yr <- 2017 + i
num_oat_uncalib[i] <- sum(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
num_oat_calib[i] <- sum(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
year_mon_cycle_tbl$mon == 6] + 1,
c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}
oat_target_uncalib_tbl <- calib_target_tbl %>%
select(year, target, group) %>%
rename(target_val = target,
target = group) %>%
filter(target %in% c("total_oat")) %>%
mutate(target = ifelse(target == "total_oat",
"Total OAT", NA),
group = "Target") %>%
bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
mutate(target = "Total OAT",
year = c(2018:2021),
group = "Uncalibrated Model")) %>%
bind_rows(., data.frame(target_val = num_oat_calib) %>%
mutate(target = "Total OAT",
year = c(2018:2021),
group = "Calibrated Model"))
oat_target_uncalib_tbl$group <- factor(oat_target_uncalib_tbl$group,
labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
levels = c("Target", "Calibrated Model", "Uncalibrated Model"))
p5 <- ggplot() +
geom_point(data = oat_target_uncalib_tbl,
aes(x = year, y = target_val, color = factor(group),
fill = factor(group))) +
xlab("Year") + ylab("Total OAT counts") +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "")
# +
# facet_wrap(~target, scales = "free") + theme(legend.position="bottom")
plotly::ggplotly(p5)Code
cbind.data.frame(parameter = calib_param_tbl_ntd$var_params,
basevalue = calib_param_tbl_ntd$basevalue,
opt_map = opt_params_notimedepend$par) %>%
filter(basevalue != opt_map) parameter basevalue opt_map
1 p__BN_PN__BPO_MISUSE 0.00032511 0.00220011
2 p__BPO_CHRONIC__BO_OD_RX 0.00097898 0.00160398
3 p__BPO_PALLIATIVE__BO_DEATH 0.37500000 0.38250000
4 p__BS_OAT_MAINT__BS_OAT_MAINT 0.70000000 0.77000000